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Abstract — This  paper  presents  the  qualitative  nature  of  com¬ 
munication  network  operations  as  abstraction  of  typical  ther¬ 
modynamic  parameters  (e.g.,  order  parameter,  temperature,  and 
pressure).  Specifically,  statistical  mechanics-inspired  models  of 
critical  phenomena  (e.g.,  phase  transitions  and  size  scaling)  for 
heterogeneous  packet  transmission  are  developed  in  terms  of  mul¬ 
tiple  intensive  parameters,  namely,  the  external  packet  load  on  the 
network  system  and  the  packet  transmission  probabilities  of  het¬ 
erogeneous  packet  types.  Network  phase  diagrams  are  constructed 
based  on  these  traffic  parameters,  and  decision  and  control  strate¬ 
gies  are  formulated  for  heterogeneous  packet  transmission  in  the 
network  system.  In  this  context,  decision  functions  and  control 
objectives  are  derived  in  closed  forms,  and  the  pertinent  results  of 
test  and  validation  on  a  simulated  network  system  are  presented. 

Index  Terms — Communication  network,  heterogeneous  packet 
transmission,  phase  transition,  statistical  mechanics. 

I.  INTRODUCTION 

CONCEPTS  of  statistical  mechanics  have  been  extensively 
used  to  model  thermodynamic  characteristics  of  physical 
phenomena  in  terms  of  their  relationships  between  micro- 
and  macrobehaviors  [1],  [2].  From  this  perspective,  tools  of 
statistical  mechanics  are  appropriate  for  modeling  the  behavior 
of  multiagent  systems  (e.g.,  communication  networks)  because 
their  global  behavior  emerges  from  the  local  dynamics  of 
the  participating  agents.  This  phenomenon  has  inspired  many 
researchers  to  investigate  complex  networks  from  the  mathe¬ 
matical  perspectives  of  thermodynamics  and  statistical  mechan¬ 
ics  [3].  Specifically,  Hui  and  Karasan  [4]  have  addressed  the 
thermodynamic  formalism  of  communication  networks.  More 
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recently,  Albert  and  Barabasi  [5]  have  reported  a  comprehen¬ 
sive  review  of  the  recent  literature  in  this  field. 

The  behavior  and  topological  organization  of  communication 
networks  [6]  and  sensor  networks  [7]  have  similar  characteris¬ 
tics  as  many  natural  and  human-engineered  systems,  such  as 
those  found  in  the  disciplines  of  sociology  [8],  biology  [9],  and 
finance  [10].  Phase  transition  is  a  characteristic  phenomenon  of 
complex  systems,  consisting  of  interacting  and  interdependent 
dynamics,  where  a  nonsmooth  change  in  the  output  behavior 
may  take  place  with  a  relatively  small  variation  of  the  system 
parameter(s).  In  this  context,  tools  of  statistical  mechanics  are 
useful  for  characterization  of  critical  phenomena  corresponding 
to  the  dependence  of  the  global  behavior  (e.g.,  connectivity, 
average  rate  of  change  of  queue  length,  and  average  packet 
drop  rate)  of  communication  networks  on  their  local  parameters 
(e.g.,  communication  radius,  packet  load,  and  transmission 
probability)  [11]-[14]. 

A  key  task  in  the  analysis  of  phase  transitions  is  the  identifi¬ 
cation  of  the  system  behavior  in  the  vicinity  of  a  critical  point, 
where  a  global  parameter  quantifies  the  presence  of  order  in 
the  underlying  system.  Usually,  this  order  parameter  [1],  [2] 
is  zero  in  the  disordered  phase  and  may  have  increased  to  a 
significant  nonzero  value  in  the  ordered  phase.  In  essence,  a 
phase  transition  is  realized  as  a  discontinuity  in  the  zeroth  or 
a  higher  derivative  of  the  order  parameter  for  a  change  from 
zero  to  a  nonzero  value  if  an  intensive  parameter  (e.g.,  tempera¬ 
ture  T)  of  the  system  is  perturbed  from  its  critical  value.  In  gen¬ 
eral,  the  nature  of  phase  transition  is  classified  into  two  types: 
1)  first-order  phase  transition,  where  the  order  parameter 
changes  discontinuously  with  the  intensive  parameter  at  the 
critical  point  and  2)  continuous  (also  called  second  or  higher  or¬ 
der)  phase  transition,  where  the  order  parameter  varies  contin¬ 
uously  with  an  intensive  parameter  during  the  phase  transition 
but  its  first  or  a  higher  derivative  with  respect  to  the  intensive 
parameter  at  the  critical  value  is  discontinuous  [15].  For  exam¬ 
ple,  ferromagnetism  is  a  well-known  case  of  continuous  phase 
transition,  where  the  order  parameter  is  magnetization  M(T ) 
with  a  nonzero  value  leading  to  spontaneous  magnetization  at 
temperatures  below  the  critical  temperature,  called  the  Curie 
point  [1],  [2]. 

The  congestion  phenomenon  in  communication  networks  is 
an  example  of  a  phase  transition.  Typically,  the  concept  of 
bottleneck  buffers  (routers)  is  used  to  detect  and  control  (i.e., 
mitigate  or  prevent)  congestion  in  communication  networks. 
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The  network  structure  could  be  reduced  to  a  multisource  and 
single  destination  (or  vice  versa )  through  the  usage  of  bottle¬ 
neck  buffers,  and  thus,  mean-field  models  have  been  developed 
for  such  cases  [16].  However,  analysis  of  certain  problems 
(e.g.,  distributed  decision  making,  perimeter  surveillance  by 
static  sensor  networks,  and  statistical  estimation  of  parameter 
distribution  over  the  entire  space)  is  apparently  very  difficult 
and  often  intractable.  Moreover,  in  many  applications,  on¬ 
line  adaptation  of  the  packet  routing  sequence  is  required  to 
maintain  acceptable  performance  in  the  presence  of  uncertain 
exogenous  perturbations  and  channel  fluctuations.  For  example, 
sensor  networks  are  often  employed  in  environments  with 
limited  communication  capability.  In  such  applications,  the  in¬ 
formation  travels  between  sensor  nodes  through  a  combination 
of  relay  nodes  and  other  sensor  nodes.  For  communication- 
constrained  sensor  networks,  the  information  that  is  transmitted 
is  usually  event  driven  (e.g.,  detection  of  intruders)  and  is 
thus  intermittent  in  nature.  For  this  reason,  it  is  important  to 
have  a  controllable  routing  procedure  that  takes  advantage  of 
the  available  network  structure  to  avoid  bottlenecks  and  other 
congestion  issues. 

From  the  aforementioned  perspective,  this  paper  presents 
a  statistical  mechanics-inspired  concept,  namely,  equilibrium 
thermodynamics  formalism,  that  is  more  suitable  in  event- 
driven  situations  compared  to  the  traditional  time-driven  mod¬ 
eling  approaches.  Phenomena  of  phase  transitions  and  size 
scaling  are  characterized  in  large-scale  communication  net¬ 
works.  The  underlying  theories  have  been  tested  and  validated 
on  a  simulation  test  bed  that  consists  of  a  2-D  square-grid 
network  model.  The  objective  here  is  to  obtain  an  unambiguous 
understanding  of  the  critical  behavior  on  a  relatively  simple 
network  with  tractable  dynamic  properties.  The  simulation 
experiments  would  help  establishing  a  knowledge  base  for 
the  synthesis  of  decision  and  control  laws  in  different  classes 
(e.g.,  wired  and  wireless)  of  real-life  communication  networks. 
Network  phase  diagrams  are  constructed  based  on  the  network 
parameters  to  formulate  control  strategies  for  heterogeneous 
packet  transmission.  In  this  regard,  closed-form  solutions  for 
a  simulated  network  system  are  presented. 

This  paper  is  organized  in  five  sections,  including  the  present 
one.  Section  II  describes  the  model  of  the  network  along  with 
the  architecture  of  the  simulation  test  bed.  Then,  the  qualitative 
nature  of  phase  transition  in  the  underlying  system  is  charac¬ 
terized,  and  the  effects  of  network  size  are  analyzed.  In  these 
analyses,  phase  transitions  are  investigated  based  on  a  single 
intensive  parameter,  namely,  the  external  packet  load  of  the 
communication  network  system.  Section  III  investigates  the  ef¬ 
fects  of  two  intensive  parameters  on  the  network  performance. 
Phase  diagrams  are  constructed  for  these  cases  by  defining 
network  analogs  of  thermodynamic  quantities  (e.g.,  order  pa¬ 
rameter,  temperature,  and  pressure).  Based  on  the  equilibrium 
phase  diagrams,  Section  IV  presents  a  novel  thermodynamic 
approach  for  controlling  heterogeneous  packet  transmission.  It 
also  presents  the  basic  philosophy  of  the  control  approach  along 
with  a  representative  numerical  experiment  on  the  network 
architecture  of  the  current  study.  The  paper  is  summarized  and 
concluded  in  Section  V  along  with  recommendations  for  future 
research. 


Fig.  1.  Network  structure  in  the  simulation  model. 


II.  Network  Simulation  Test  Bed 

The  network  test  bed  simulates  a  2-D  square  grid  [12]  as 
shown  in  Fig.  1,  where  the  nodes  (routers)  are  placed  at  the 
grid  points.  For  a  square-grid  network  with  TV  x  N  nodes,  there 
are  (47V  —  4)  boundary  nodes  (shown  as  squares  in  Fig.  1)  and 
(TV2  —  4 AT  +  4)  internal  nodes  (shown  as  circles  in  Fig.  1). 
Only  boundary  nodes  are  assumed  to  be  the  sources  and/or 
the  sinks  for  packet  generation  and  destruction;  internal  nodes 
can  only  transmit  the  packets.  Each  node  receives  packets  in  a 
queue  (with  a  finite  buffer  length)  from  its  neighboring  nodes, 
and  packets  are  terminated  after  reaching  their  respective  desti¬ 
nations.  In  each  time  unit,  packets  are  created  in  the  boundary 
nodes  with  a  Poisson  arrival  rate  A.  For  simplicity,  all  packets 
are  considered  to  be  of  unit  length.  The  destination  of  a  packet 
is  chosen  randomly  from  the  set  of  boundary  nodes,  including 
its  source  node.  Each  node  (boundary  or  internal)  transmits 
one  packet  from  the  head  of  its  queue  to  a  deterministically 
chosen  neighboring  node  at  each  time  unit.  The  node  chosen 
to  forward  a  data  packet  is  selected  so  that  the  packet  travels 
via  the  shortest  path  to  its  destination.  When  there  are  more 
than  one  candidate  nodes  for  the  shortest  path,  the  node  with  a 
smaller  queue  length  is  chosen  to  reduce  the  probability  of  early 
congestion  in  the  network.  A  node  drops  the  oldest  packet  from 
its  fully  occupied  queue  to  accommodate  a  new  packet. 


A.  Network  Congestion  Modeled  as  a  Phase  Transition 

Congestion  is  a  phenomenon  of  significant  importance  in 
communication  networks.  Before  a  congestion  occurs  (i.e.,  the 
network  is  capable  of  handling  the  average  number  of  arriving 
packets),  the  average  packet  drop  rate  remains  negligibly  small. 
However,  as  the  packet  influx  rate  crosses  a  critical  threshold, 
the  drop  rate  increases  to  a  nonnegligible  value.  Thus,  at  a 
steady  state,  the  nonzero  average  packet  drop  rate  can  serve 
as  an  index  of  the  degree  of  congestion.  The  phenomenon  of 
congestion  can  be  viewed  as  a  continuous  phase  transition  from 
a  steady  state  to  an  unsteady  state  of  network  communications, 
which  is  characterized  based  on  the  concepts  of  equilibrium 
thermodynamics  in  this  paper.  Previous  papers  [11],  [17]  re¬ 
ported  similar  characterization  based  on  the  average  rate  of 
change  of  queue  length;  however,  they  required  the  assumption 


SARKAR  etal.:  MODELING  OF  HETEROGENEOUS  PACKET  TRANSMISSION 


1085 


of  infinite  queue  length.  This  assumption  is  removed  in  this 
paper  to  address  the  congestion  problem  in  a  more  realistic 
scenario. 

As  discussed  in  Section  I,  a  global  order  parameter  of  the  sys¬ 
tem  under  consideration  needs  to  be  identified  for  investigating 
the  occurrence  of  phase  transitions  in  communication  networks. 
The  packet  drop  rate  is  a  feasible  candidate  to  serve  as  the 
order  parameter  in  this  analysis.  Furthermore,  global  intensive 
parameters  that  trigger  the  network  phase  transition  need  to  be 
identified,  and  in  this  problem,  packet  influx  rate  A  is  one  such 
parameter.  Since  only  4 (TV  —  1)  boundary  nodes  generate  the 
packets  in  this  simulated  network  and  all  N 2  nodes  receive 
these  packets,  a  surface  correction  is  needed  to  accurately 
define  the  intensive  parameter,  namely,  effective  load  per  node 
A eff,  which  is  defined  as  follows  to  avoid  the  introduction  of 
the  surface  effects: 

A  few  scaling  adjustments  are  made  to  define  the  order 
parameter.  Let  d(t)  be  the  total  number  of  packets  dropped  in 
the  network  at  time  t.  Thus,  the  instantaneous  packet  drop  rate 
(per  node)  at  time  t  is  d(t)/N 2  and  its  time-averaged  value  is 

denoted  as  D  =  d(t)/N2.  The  order  parameter  M  is  defined 
based  on  the  normalization  of  D  with  respect  to  A ef  / 

M=TL.  (2) 

Xeff 

Note  that  the  order  parameter  M  is  the  fraction  of  incoming 
packets  that  are  dropped  from  the  network  nodes  per  unit  time. 
Therefore,  in  the  absence  of  congestion,  M  could  be  assumed  to 
be  negligibly  small,  i.e.,  there  is  no  packet  drop.  As  the  network 
becomes  congested,  the  worst  case  scenario  is  the  dropping  of 
all  incoming  packets,  i.e.,  M  approaches  1. 

Following  the  aforementioned  procedure,  M  is  computed  for 
given  values  of  A eff-  To  eliminate  the  effects  of  transients, 
the  expected  value  of  D  is  calculated  using  only  steady- state 
time  series  data.  The  plot  of  M  versus  Ae//  for  the  10  x  10 
simulated  network  is  shown  in  Fig.  2.  It  is  seen  that  there  is  a 
critical  value  Xce ^  «  0.11  of  effective  load  per  node  such  that, 
for  A  eff  <  A  Iff,  the  order  parameter  M  is  almost  negligible; 
in  contrast,  for  A eff  >  A M  takes  on  nonzero  values.  This 
change  of  network  behavior  across  the  critical  value  Xeff  1S 
identified  as  a  continuous  phase  transition ,  where  the  network 
moves  from  an  uncongested  phase  (U)  of  negligibly  small  M 
to  a  congested  phase  (C)  of  finite  positive  M. 

B.  Construction  of  Size  Scaling  Laws 

It  is  very  important  to  understand  the  laws  of  size  scaling  for 
human-engineered  multiagent  systems  (e.g.,  communication 
network  systems)  due  to  their  inherent  finite  and  smaller  sizes, 
as  compared  to  their  natural  counterpart.  In  general,  for  finite¬ 
sized  systems,  the  critical  value  of  the  intensive  parameter 
Agjj(TV)  is  a  function  of  size  N,  and  as  N  goes  to  infin¬ 
ity,  Agjj(TV)  converges  to  the  Agjj(oo)  of  the  corresponding 


Fig.  2.  Continuous  phase  transition  in  a  square-grid  communication  network. 
The  normalized  packet  drop  rate  (M)  is  plotted  against  the  effective  load  per 
node  (Ae//). 

infinite- sized  system  in  the  thermodynamic  limit.  In  the  sta¬ 
tistical  mechanics  literature,  the  space  correlation  length  £(A) 
behaves  as  a  function  of  |A|  in  a  power  law.  Here,  A  can  be 
considered  as  A ce^{N)  —  Xce^{oo).  Also,  the  following  size 
dependence  form  of  a  test  function  £(A)  is  assumed  [18],  which 
conforms  to  the  characteristics  of  £,  and: 

t(\ceff(N)-\ceff(oo))=PN  (3) 

where  p  is  the  proportionality  constant.  Thus,  the  following  law 
is  postulated  based  on  the  structure  of  the  2-D  Ising  model  (also 
known  as  the  Onsager  model)  [1],  [2]: 

Kff(N)  =  X^ff(oo)  +  qN~”  (4) 

where  q  and  v  are  constant  model  parameters  that  are  called 

the  coefficient  and  index,  respectively.  However,  in  the  present 
problem,  Ae//  ^  (1/7V).  Thus,  it  is  impossible  to  provide  a 
finite  effective  load  per  node  to  an  infinite- sized  network,  which 
implies  that  there  cannot  be  any  finite  load  phase  transition  in  an 
infinite- sized  network,  i.e.,  Ag^(oo)  =  0.  Hence,  for  this  case 

Xceff(N)=qN-i  (5) 

where  the  model  parameters  (i.e.,  coefficient  q  and  exponent  v) 
in  (5)  are  identified  by  fitting  with  the  simulation  data  generated 
on  the  network  test  bed.  Fig.  3  shows  the  closeness  of  the  fitted 
data  with  the  model,  where  the  identified  parameters  are  as 
follows:  q  =  100  034  «  1.08  and  v  =  (1/1.1)  «  0.9. 

Fig.  4  shows  the  possibility  of  generating  a  size-independent 
global  model  for  phase  transition  in  square-grid  communication 
networks.  Such  a  model  can  be  constructed  by  using  a  reduced 
(i.e.,  normalized)  effective  load  A red  =  \eff/\ce ^  instead  of 
Xeff .  This  approach  draws  inspiration  from  classical  ther¬ 
modynamics,  where  compressibility  curves  for  different  pure 
substances  are  unified  by  using  reduced  pressure  or  reduced 
temperature  instead  of  directly  using  pressure  or  temperature 
as  the  independent  variables. 
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Fig.  3 .  Size  dependence  of  critical  network  load.  The  critical  value  of  effective 
load  per  node  (A£ ^ )  is  plotted  against  the  size  of  the  network  ( N )  in  the  log- 
log  scale. 


Fig.  5.  Average  (over  time)  packet  drop  rate  (D)  as  a  function  of  packet 
arrival  rate  (A eff)  and  transmission  probability  (P). 

III.  Heterogeneous  Packet  Transmission: 
Introduction  of  a  Second  Intensive  Parameter 


Fig.  4.  Size-independent  global  model  for  network  phase  transition.  The 
normalized  packet  drop  rate  (M)  is  plotted  against  the  effective  load  per  node 
(A eff),  normalized  with  respect  to  its  critical  value  (A 

Remark  2.1  ( Critical  Slowing  Down):  Note  that  there  is  a 
slight  mismatch  among  the  phase  transition  curves  in  Fig.  4, 
which  can  be  attributed  to  the  fact  of  critical  slowing  down  [18]. 
The  normalized  time-correlation  function  (j>x(ti  —  £2)  of  an 
observable  x(t)  is  defined  [19]  as 


Previous  analyses  studied  the  effects  of  variations  of  a  single 
parameter  (e.g.,  packet  arrival  rate  A)  on  the  average  packet 
drop  rate  in  the  network.  This  section  introduces  a  second 
network  parameter,  namely,  the  packet  transmission  probability 
P,  for  networks  transmitting  heterogeneous  packets.  Given  that 
each  node  of  the  network  has  multiple  independent  queues 
(buffers)  for  different  packet  types,  packet  transmission  prob¬ 
ability  is  defined  as  the  probability  of  transmitting  a  particular 
packet  in  a  time  unit. 

In  the  previous  case,  with  a  single  packet  type,  P  was  set 
to  unity.  It  is  obvious  that,  when  P  is  decreased  from  unity, 
the  network  is  expected  to  move  from  the  steady  phase  to  an 
unsteady  phase  for  lower  values  of  the  effective  load  \eff- 
Similar  phenomena  take  place  in  two-phase  thermodynamic 
(e.g.,  solid-liquid  and  gas-liquid)  systems,  where  the  critical 
temperature  of  solid-liquid  (or  gas-liquid)  transition  can  be 
altered  by  changing  the  superincumbent  pressure;  lower  pres¬ 
sures  usually  lead  to  lower  values  of  melting  or  boiling  point. 
In  this  context,  the  packet  transmission  probability  P  of  the 
nodes  is  called  the  network  pressure ,  while  A ef  f  is  called  the 
network  temperature.  Therefore,  phase  transition  in  the  network 
is  a  function  of  a  combination  of  network  temperature  Ae//  and 
network  pressure  P.  Fig.  5  presents  a  3-D  plot  of  the  average 
packet  drop  rate  as  a  function  of  A e/  /  and  P,  and  it  is  observed 
that  Agjy  decreases  with  a  decrease  in  P. 


4>x(fi  —  £2) 


[(x(£i)x(£2))  -  ( x )2]  i*i-t2i 

[<**>  -  (*)2]  *C 


(6) 


where  r  is  the  slowest  relaxation  time  or  temporal  correlation 
length  of  the  observable  x(t).  It  is  known  that,  near  a  phase 
transition  point,  the  relaxation  time  r  of  the  slowest  mode  of 
a  system  diverges,  i.e.,  r  — >  00.  This  phenomenon  is  known  as 
the  critical  slowing  down ,  and  as  a  consequence,  it  takes  a  long 
time  to  make  two  consecutive  independent  observations  near 
a  phase  transition  point.  Therefore,  it  is  difficult  to  simulate  a 
large  system  in  the  vicinity  of  the  phase  transition  point. 


A.  Problem  Formulation 

The  problem  of  heterogeneous  packet  transmission  in  net¬ 
works  is  formalized  in  this  section.  Let  there  be  K  independent 
queues  for  K  types  of  packets  for  each  node  of  the  network. 
For  example,  in  the  context  of  perimeter  surveillance  sensor 
networks,  heterogeneity  may  accrue  from  differences  in  the 
sensor  modality  (e.g.,  audio,  video,  and  magnetic  sensor  pack¬ 
ets).  In  terms  of  packet  arrival,  let  the  effective  Poisson  arrival 
rate  of  packet  type  i  be  simply  (subscript  eff  is  omitted) 
denoted  by  A^  (0  <  A^  <  1)  for  i  =  1,  2, . . .  K.  However,  only 
one  channel  is  used  for  packet  transmission,  i.e.,  at  a  given 
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Node 


Packet  Transmission 
Channel 

Fig.  6.  Illustration  for  heterogeneous  packet  transmission  mechanism. 

time  instant,  a  node  can  at  most  transmit  one  packet.  A  typical 
scenario  is  illustrated  in  Fig.  6.  The  probability  that  a  node  in 
the  network  transmits  a  packet  of  type  i  is  denoted  by  Pi,  where 
0  <  Pi  <  1.  A  single  transmission  channel  leads  to  the  con¬ 
straint  Pi  <  1-  In  this  setting,  the  network  is  analogous 
to  multicomponent  materials  in  a  thermodynamic  sense.  Thus, 
phase  diagrams  can  be  constructed  for  the  equilibrium  states  for 
such  networks  which  eventually  lead  to  a  probabilistic  (static) 
control  strategy  of  heterogeneous  packet  transmission.  The  next 
section  presents  the  idea  of  phase  diagrams  for  networks  trans¬ 
mitting  heterogeneous  packets  with  a  representative  example  of 
two  packet  types. 


Independent  Queues  for 
Heterogeneous  Packets 


B.  Construction  of  Phase  Diagrams 

In  the  current  setting  of  independent  queues  for  different 
packet  types,  the  uncongested  phase  U  (respectively,  congested 
phase  C)  for  a  particular  packet  type  is  characterized  by  zero 
(respectively,  nonzero)  packet  drop  rate.  Thus,  with  respect  to 
a  particular  packet  type,  a  network  is  in  an  uncongested  or  a 
congested  phase.  This  leads  to  the  existence  of  mixed  phases, 
where  the  network  remains  in  an  uncongested  phase  for  some 
packet  types  while  being  in  a  congested  phase  for  the  rest  of  the 
packet  types.  For  a  network  carrying  K  packet  types,  there  are 
2k  possible  phases,  which  include  the  fully  uncongested,  fully 
congested,  and  mixed  phases. 

The  purpose  of  a  network  phase  diagram  is  to  determine  its 
phase  for  given  distributions  of  A i  and  Pi.  A  representative 
example  with  two  packet  types  is  presented  here.  As  discussed 
earlier,  there  can  be  the  following  four  different  network  phases 
for  this  example. 

1)  Completely  uncongested  phase  (U1  +  U2):  Both  packet 
types  are  in  uncongested  phase,  which  signifies  negligibly 
small  values  of  average  packet  drop  rates  Di  and  D2  for 
both  queues. 

2)  Congested  packet  type  1  and  uncongested  packet  type  2 
(Cl  +  U2):  Packet  type  1  has  a  relatively  large  (i.e., 
nonzero)  drop  rate  D\  while  packet  type  2  still  maintains 
a  negligibly  small  D2 . 

3)  Uncongested  packet  type  I  and  congested  packet  type  II 
(U1  +  C2):  Packet  type  1  maintains  a  negligibly  small 
drop  rate  D\  while  packet  type  2  has  a  relatively  large 
(positive)  drop  rate  D2 . 

4)  Completely  congested  phase  (Cl  +  C2):  Both  packet 
types  are  in  the  congested  phase,  i.e.,  average  packet  drop 
rates  Di  and  D2  are  large  for  both  queues. 

Remark  3.1 :  When  compared  to  a  multicomponent  material 
in  the  thermodynamic  sense,  the  network  with  two  packet  types 
has  a  remarkable  similarity  with  two-component  mixtures,  e.g., 


CL 


U^U2 

Cl  +  U2 


Ui  +  C2 
Ci  +  C2 


Fig.  7.  Network  phase  diagram  with  two  packet  types.  The  transmission 
probability  of  packet  type  i  is  denoted  as  Pi  under  the  constraint  Pi  +  P2  =  1. 
The  arrival  rate  for  packets  of  type  i  is  denoted  at  A i . 

a  mixture  of  olivine  (i.e.,  an  isolated  tetrahedra)  and  pyroxene 
(i.e.,  single  chain  tetrahedra).  The  corresponding  phase  diagram 
is  known  as  the  binary  eutectic  phase  diagram  [20]  that  ex¬ 
plains  the  chemical  process  of  generating  the  two  immiscible 
crystals  from  a  completely  miscible  liquid  based  on  temperature 
and  pressure  variations.  This  process  also  has  four  possible 
phases,  which  are  the  following:  completely  liquid  phase  (anal¬ 
ogous  to  Cl  +  C2),  mixture  of  liquid  olivine  with  pyroxene 
crystals  (analogous  to  Cl  +  U2),  mixture  of  olivine  crystals 
with  liquid  pyroxene  (analogous  to  U1  +  C2),  and,  finally,  the 
completely  crystallized  phase  (analogous  to  U1  +  U2). 

Fig.  7  shows  a  3-D  phase  diagram  for  the  aforementioned 
example,  which  is  constructed  via  Monte  Carlo  simulation. 
Four  different  phases  are  color  coded.  Both  Ai  and  A2  are  varied 
independently  from  0  to  1.  Under  the  constraint  Pi  +  P2  =  1, 
planes  for  three  different  values  of  Pi,  namely,  0.2,  0.5,  and  0.8, 
show  changes  in  the  sizes  of  the  phases  across  the  range  of  Pi . 
The  regularity  in  shapes  of  different  phase  zones  is  attributed 
to  the  independence  of  the  queues  in  different  packet  types. 
This  type  of  phase  diagrams  leads  to  the  concept  of  controlling 
heterogeneous  packet  transmission  from  the  perspectives  of 
statistical  mechanics.  The  basic  idea  is  to  move  as  close  as 
possible  to  the  phase  boundary  of  the  U1  +  U2  phase  from 
outside  by  choosing  an  appropriate  Pi.  On  the  other  hand,  if 
the  network  is  already  inside  the  U1  +  U2  phase,  then  the 
strategy  should  be  to  move  as  far  away  as  possible  from  its 
phase  boundary,  also  by  choosing  an  appropriate  Pi.  The  key 
objective  is  to  achieve  robustness  against  the  uncertainties  in  Ai 
and  A 2  or  any  other  internal  disturbances. 


IV.  Control  of  Packet  Transmission 

The  thermodynamic  approach,  presented  earlier  to  analyze 
and  understand  the  congestion  phenomenon  in  communication 
networks,  could  lead  to  a  novel  means  to  control  heteroge¬ 
neous  packet  transmission  by  tuning  the  packet  transmission 
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probability  at  a  node  level.  This  section  presents  the  basic 
philosophy  of  such  a  control  strategy  along  with  representative 
numerical  experiments  on  the  network  test  bed.  From  a  control 
perspective,  the  congestion  states  of  a  network  are  categorized 
into  the  following  two  types. 

1)  Phase  Type  I:  At  least  a  queue  for  one  packet  type  is  in 
the  congested  phase. 

2)  Phase  Type  II:  Queues  for  all  packet  types  are  not  in  the 
congested  phase. 

The  control  strategies  are  fundamentally  different  for  these 
two  phase  types.  For  Phase  Type  I,  the  drop  rate  of  at  least  one 
type  of  packet  is  nonzero,  where  the  objective  is  to  choose  a 
probability  distribution  of  packet  transmission  to  minimize  the 
worst  case  drop  rate  among  packet  types  with  nonzero  drop 
rates  that  can  be  scaled  by  respective  packet  arrival  rates  and 
also  by  a  user-defined  packet  priority  distribution.  In  contrast, 
in  Phase  Type  II,  drop  rates  for  all  packet  types  are  already 
zero.  Thus,  any  transmission  probability  distribution  that  keeps 
all  the  drop  rates  at  zero  is  an  optimal  distribution  according 
to  the  objective  of  Phase  Type  I.  However,  a  unique  solution 
can  be  achieved  by  the  use  of  a  different  objective  that  would 
provide  robustness  against  the  uncertainties  and/or  fluctuations 
in  packet  arrival  rates.  The  idea  here  is  to  choose  a  transmission 
probability  distribution  that  pushes  the  network  deeper  into  the 
steady  phase  (i.e.,  away  from  the  phase  boundary)  so  that,  even 
with  fluctuations  or  uncertainties  in  the  packet  arrival  rate,  the 
probability  of  being  in  the  unsteady  phase  is  minimized.  In 
other  words,  the  objective  is  to  maximize  the  least  uncertainty 
tolerance  in  packet  arrival  rate.  As  before,  the  uncertainty 
tolerance  can  be  scaled  by  the  respective  packet  arrival  rates 
and  the  user-defined  packet  priority  distribution.  In  the  sequel, 
the  control  objectives  are  described  formally  along  with  optimal 
solutions  for  general  and  special  cases. 


The  optimal  packet  transmission  probability  for  each  packet 
type  i  is  denoted  as  P*  and  is  given  by  the  constrained  opti¬ 
mization  solution  as 


(Pi,  AV  •  •  P*k)  =  arg  min  C(PU  P2,...PK) 

E,  ^ 

<Xifl(Pi,\i) 

=  arg  mm  max - - - . 

V.P4<1  i 


(10) 

(11) 


Although  the  constraint  on  the  transmission  probability  dis¬ 
tribution  is  JA  Pi  <  1,  the  inequality  is  replaced  by  an  equal¬ 
ity  sign  for  optimality.  The  rationale  for  involving  Xz  in  the 
objective  function  is  to  normalize  the  drop  rates  for  each 
packet  type  as  it  was  done  to  define  the  order  parameter  M 
in  Section  II.  In  this  context,  the  solution  of  the  aforementioned 
optimization  problem  is  fairly  straightforward.  For  the  optimal 
packet  transmission  probability  distribution,  a  solution  of  the 
form 


OiiDi  djDj 


Vi,  j 


(12) 


is  globally  optimal  if  it  exists.  This  is  because,  if  one  chooses  a 
packet  transmission  probability  distribution  to  reduce  a^Dij \ 
even  further  for  some  value  of  i ,  then  ajDj/Xj  would  increase 
for  at  least  one  j.  Thus,  the  cost  functional  would  increase. 
The  proof  of  optimality  is  presented  in  Appendix  A  which 
also  introduces  a  sequential  and  recursive  approach  to  min-max 
optimization  when  the  condition  described  in  (12)  does  not 
exist. 

2 )  Control  Objective  for  Phase  Type  II:  In  Phase  Type  II, 
Di  =  0,  i.e.,  /i (Pi,  Xi)  =  0  Vi,  which  leads  to  the  critical  curve 
that  yields  the  critical  transmission  probability  Pc  for  a  given 
packet  arrival  rate  A.  Let  the  critical  curve  be  described  as 


A.  Control  Objectives  and  Optimal  Solutions 

A  user-defined  packet  priority  distribution  {c^},  with  oti  G 
[0, 1]  and  ai  =  is  considered  for  defining  the  control 
objectives  in  both  phase  types.  A  higher  value  of  oti  signifies 
higher  priority  for  packet  type  i  while  lowering  the  drop  rate  in 
Phase  Type  I  or  increasing  the  uncertainty  tolerance  in  Phase 
Type  II. 

I )  Control  Objective  for  Phase  Type  I:  In  Phase  Type  I,  the 
following  model  for  the  average  drop  rate  Di  of  packet  type  i  is 
assumed: 


Di  =  f1(Pi,Xi)  Vi  (7) 


where  Pi  and  A*  are  defined  as  before.  To  find  an  optimal 
transmission  probability  distribution  Pi  for  a  given  packet 
arrival  distribution  Xt  and  a  packet  priority  distribution  cu  that 
minimize  the  (weighted)  worst  case  packet  drop  rate,  the  cost 
functional  is  chosen  as 


C(P1,P2,...PK) 


a,D, 
max  — - — 
i=l,...,K  A  i 


:/l(A,V) 
A 


(8) 

(9) 


Pc  =  /2(  A).  (13) 

For  a  given  packet  arrival  distribution  A*,  if  the  network  is  in 
Phase  Type  II,  then 

^ ^  L (A, •)<  1.  (14) 

i  i 

This  implies  that  the  packet  transmission  probability  distri¬ 
bution  Pf,  where  Pj  <  1),  is  sufficient  for  the  network 
to  be  in  Phase  Type  II.  However,  the  queues  for  different 
packet  types  can  still  use  the  remaining  packet  transmission 
capability  of  the  network  expressed  as  (1  —  JA  Pf).  The  goal 
here  is  to  distribute  this  remaining  packet  transmission  capabil¬ 
ity  among  the  queues  in  an  optimal  sense.  To  formalize  this 
problem,  a  virtual  packet  arrival  distribution  Xf  is  defined, 
where  Xf  (Xi  <  Xf  <  1)  denotes  a  distribution  that  includes 
uncertainty/fluctuations  in  packet  arrival  rates,  for  which  the 
network  just  remains  in  Phase  Type  II.  The  corresponding 
packet  transmission  probability  Pi  is  given  as  f2(Xf),  which 
leads  to 


Af  —  A^  +  A^  (15) 

=^A  i=/2-1(Pi)-Ai  (16) 
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where  Ai  (>  0)  uncertainty  tolerance  provided  for  packet  type  i, 
for  which  the  network  still  remains  in  the  steady  phase  (Phase 
Type  II).  The  objective  of  the  optimization  problem  for  Phase 
Type  II  is  to  choose  a  transmission  probability  distribution  Pi 
that  maximizes  the  minimum  uncertainty  tolerance  for  packet 
arrival  statistics.  Thus,  the  objective  functional  is  the  worst  case 
(weighted)  tolerance  Ai  that  is  given  as 


0(P1,P2,...,PK) 


.  A  i 

mm  — 

OLi 


min 

i 


f£\PQ-Xi 


OLi 


(17) 

(18) 


The  constrained  optimization  problem  that  has  to  be  solved 
to  obtain  the  optimal  packet  transmission  probability  is  given  as 


(P^,Pi,...P*K)  =  arg  max  0(P1,P2, . . .  PK)  (19) 

y.PjKi 

Pi>f2W 

=  arg  max  mm  — - .  (20) 

EiP*^  *  ai 

Pi>f2W 

The  solution  for  the  aforementioned  optimization  problem  is 


A  i  Xj 

ai  aj 


Vi,  j. 


(21) 


The  details  are  presented  in  Appendix  B.  Note  the  use  of  the 
packet  priority  distribution  ai  in  the  optimization  problem, 
where  the  queue  for  packet  type  i  gains  more  uncertainty 
tolerance  for  a  higher  value  of  ai. 


B.  Approximation  of  Functional  Forms 

The  formulation  of  a  general  control  strategy  is  presented 
above.  For  a  given  network  structure  and  routing  strategy,  one 
needs  to  estimate  the  functions  /i  and  /2,  for  Phase  Type  I  and 
Phase  Type  2,  respectively,  in  order  to  implement  the  control 
strategy  (see  Section  IV- A).  An  example  procedure  is  shown 
in  this  section  with  respect  to  the  network  considered  in  the 
current  study.  Fig.  5  generates  an  idea  for  constructing  the 
functional  forms  for  both  /i  and  For  Phase  Type  I  (i.e., 
average  packet  drop  rate  per  node  D  >  0),  the  function  /i  can 
be  approximated  with  a  linear  function,  i.e.,  with  a  plane.  As 
there  is  no  bias,  the  plane  equation  for  D  as  a  function  of  P 
and  A  is 


D  =  aP  +  bX.  (22) 

The  linear  least  square  fit  of  the  (numerical)  experimental  data 
is  shown  in  Fig.  8,  and  the  error  is  found  to  be  sufficiently  low. 
Note  that  the  plane  (22)  only  fits  the  data  in  Phase  Type  I  and  is 
not  valid  for  Phase  Type  II  data,  where  D  is  always  zero.  The 
error  in  the  fit  is  observed  to  increase  near  the  critical  curve 
(see  the  circled  zone  in  the  figure).  This  can  be  attributed  to 
the  problem  critical  slowing  down  near  the  critical  points  as 
discussed  earlier.  Note  that  D  decreases  with  an  increase  in  P 
and  D  increases  with  an  increase  in  A. 


Q 


Fig.  8.  Approximation  of  average  (over  time)  packet  drop  rate  ( D )  as  a 
function  of  packet  arrival  rate  (A e//)  and  transmission  probability  (P)  in 
square-grid  networks.  (Note:  The  region  for  D  <  0  is  inaccessible  because  the 
packet  drop  rate  is  nonnegative.) 


The  critical  curve  equation  (for  Phase  Type  II)  is  found  by 
setting  D  =  0.  In  other  words,  the  critical  curve  is  the  straight 
line  at  which  the  least  square  fit  plane  intersects  the  D  =  0 
plane.  Thus,  the  following  equation  for  the  critical  curve  is 
obtained: 


Pc  =  cX 


(23) 


where  c  is  a  positive  constant  as  Pc  must  increase  with  an 
increase  in  A.  Furthermore,  c  is  expressed  as  follows:  c  = 
—b/a,  where  b  is  positive  and  a  is  strictly  negative  (i.e.,  a  = 
—  \a\  <  0).  Fig.  8  shows  the  critical  line  (dotted)  for  the  current 
network  model. 

The  aforementioned  approximate  functional  forms  allow  the 
derivation  of  the  closed-form  solutions  (for  both  phase  types) 
of  the  optimal  packet  transmission  probability  distribution  by 
solving  algebraic  equations. 

1 )  Solution  for  Phase  Type  I:  In  Phase  Type  I,  the  suffi¬ 
ciency  criterion  for  the  optimal  distribution  is  given  by  the 
solution  of  the  following  system  of  algebraic  equations: 

—Pi+aib=y—Pj+ajb  \/i,  j  G  (1, . . . ,  K}.  (24) 

A  i  Aj 


For  a  network  with  K  packet  types,  the  above  system  pro¬ 
vides  (K  —  1)  independent  equations,  and  the  K th  equation  is 
provided  by  the  constraint  JA  P{  =  1. 

Defining  Q  =  (< c^a/A*)P^  +  aib  Vi  G  {1, ... ,  K },  (23)  and 
(24)  are  combined  to  obtain 

Q=^M(pc_pt)  Vi  G  (1, ,  K}.  (25) 

Ai 

It  follows  that  Pi  can  be  evaluated  in  terms  of  Q  as 


(26) 
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The  parameter  Q  is  evaluated  by  summing  the  expressions  for 
all  i  and  using  the  constraint  ^  •  Pi  =  1 


network  with  two  packet  types,  the  optimal  packet  transmission 
probability  for  type  1  is  given  by  Pi 


K 


EPf- 1 


(27) 


fo 


Pi- 


q2+cA2(qi-q2) 

Mm] 


if  Ol  <  max 

a2 


1-^r.O 


<  <  _ l _ 

“  “2  “  max[l-^.o] 

if  ^  — T 


ifmax 


However,  for  certain  choices  of  c**,  the  solution  of  Pi  may  be 
negative  for  a  given  distribution  of  A i,  which  does  not  satisfy 
the  constraint  Pi  >  0  Vi.  This  implies  that  a  solution  for  the 
set  of  equations  in  (24)  does  not  exist  in  the  feasible  region. 
Therefore,  as  described  in  Appendix  A,  the  negative  solutions 
should  be  constrained  to  zeros,  and  the  optimization  should 
be  repeated  for  the  rest  of  the  packet  types  in  order  to  satisfy 
the  constraint  JA  Pi  —  1.  The  algorithm  for  such  a  sequential 
optimization  procedure  is  given  hereinafter. 


and  P2  =  1  —  Pi-  Note  that,  in  Phase  Type  I,  J2i  Pi  — 
cJ2i^  > 1  which  ensures  that  max[l  —  (1/ cA2 ) ,  0]  < 
(1/ max[l  -  (l/cAi),0]). 

2 )  Solution  for  Phase  Type  II:  For  Phase  Type  II,  as  dis¬ 
cussed  earlier,  the  remaining  transmission  capability  (beyond 
the  critical  curve)  1  —  JA  P?  is  distributed  among  the  queues 
of  different  packet  types.  To  achieve  the  closed-form  solution, 
we  begin  with  the  solution  for  Phase  Type  II  [given  in  (21)]. 
Using  the  model  in  (23),  the  optimal  solution  can  be  written  as 


Algorithm  for  Sequential  Optimization 

0(1)  =  (Ef=i^-i)/Ef=i(V«fc|«l); 

Pf  =  Pf  -  (Q^ \j / ctj\a\);  (Evaluate  pf  Vj) 

i  =  1; 

while  Pj  <  0  for  some  jdo 

for  all  j  :  Pf  <  Odo 

Pf  <-  0; 

end  for 

i  i  +  1 

Q{i)  =  (Ek={j..Pr^0}P^-i)/ 

for  all  j  :  Pf~1]  ^  Odo 

/f  =J-_(QWA./a.|a|); 

end  for 

end  while 

Note  that,  as  a  consequence  of  the  aforementioned  algorithm, 
the  transmission  probability  distribution  may  have  Pi  =  0  for 
some  values  of  i  and  the  rest  of  the  elements  will  have  Pj  >  0 
(where  j  G  {1,  2, . . . ,  K}  \  {i}).  The  scaled  drop  rates  can  be 
written  as 


oaDi 


OLih  - 


OLi\a\ 


Pi • 


(28) 


=  R 


V  i 


(30) 


where  R  is  a  constant.  Appendix  B  establishes  the  existence  and 
optimality  of  the  solution.  Therefore 


Pi  —  c(Rai  +  Xi).  (31) 


The  value  of  constant  R  is  evaluated  using  the  constraint 

T,kpk  =  cT,k(Rak  +  V)  =  1  as 

1  -Va,..  (32) 

c  ^ 


Thus,  the  optimal  packet  transmission  probability  distribution 
for  Phase  Type  II  is  given  by 


Pi  —  c\i  -f  OLi 


1  -  A k 

k 


(33) 


For  example,  the  exact  solution  for  the  network  with  two  packet 
types  is 


Pi  —  +  c  [A i  —  ai(Xi  +  A2)]  (34) 


and  P?  =  1  —  P\ . 


C.  Representative  Experimentation 


Therefore,  Vj,  k  e  {1,  2, . . . ,  K}  \  {i}  (i.e.,  for  all  packet  types 
with  nonzero  transmission  probability),  the  solution  obtained 
from  the  algorithm  will  satisfy 


^3  ^  3  ^  /c  Dk 

X  j  X  k 


(29) 


Thus,  due  to  the  constraints  on  the  transmission  probability 
distribution,  the  optimal  condition  described  in  (12)  will  only 
be  satisfied  on  a  subset  (containing  packet  types  with  nonzero 
transmission  probability)  of  the  set  of  all  packet  types.  The 
rest  of  the  packet  types  will  have  zero  transmission  probability 
essentially  due  to  their  low  ( a/X )  ratio.  As  an  example,  for  a 


A  representative  experimentation  on  the  simulated  network 
test  bed  is  presented  in  this  section  to  validate  the  effectiveness 
of  the  control  strategy  developed  earlier.  The  simulation  test  bed 
is  the  same  network  system  analyzed  in  this  paper.  This  network 
structure  has  the  potential  of  serving  as  a  fundamental  building 
block  for  real-life  (and  hence  more  complex)  networks.  For 
example,  the  network  structure  used  here  can  be  considered 
as  a  part  of  a  larger  sensor  network,  where  edges  of  the 
current  (smaller)  network  interact  with  the  rest  of  the  network. 
In  [21],  an  important  observation  is  made  that  distinguishes 
sensor  network  routing  from  other  classical  IP-based  routing 
schemes.  In  sensor  networks,  generated  data  have  significant 
redundancy  as  multiple  sensors  in  the  vicinity  of  an  event  may 
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Fig.  9.  System  trajectory  in  the  phase  diagram  for  the  representative 
experimentation.  The  transmission  probability  of  packet  type  i  is  denoted  as 
Pi  under  the  constraint  Pi  +  P2  =  1.  The  arrival  rate  for  packets  of  type  i  is 
denoted  at  A  i . 

generate  the  same  data.  Routing  protocols  should  use  such 
redundancy  to  improve  energy  and  bandwidth  utilization.  The 
redundancy  can  be  used  by  implementing  an  activity  schedule 
where  a  significant  fraction  of  nodes  are  kept  inactive  at  a  given 
time.  Furthermore,  often  events  to  be  detected  by  large  sensor 
networks  are  local,  rare,  and  random.  Such  information  needs 
to  be  propagated  throughout  the  network  for  activity  scheduling 
as  well  as  for  achieving  robustness  to  node  failures.  This  paper 
addresses  the  issue  of  information  propagation  by  multisource 
and  multidestination  types  of  communication. 

Two  different  types  of  packets,  namely,  packet  type  1  and 
packet  type  2  (e.g.,  video  and  audio  sensors),  are  considered 
for  the  validation  of  the  underlying  concept.  Following  the 
notations  used  in  this  paper,  the  packet  arrival  statistics  is 
described  by  Ai  and  A2,  packet  priority  distribution  is  denoted 
by  a  1  (taken  as  0.6)  and  a2  (taken  as  0.4),  and  drop  rates  are 
denoted  by  Di  and  D2.  The  goal  here  is  to  estimate  Ai  and 
A2  over  a  time  window,  then  to  detect  the  phase  type  with 
respect  to  the  estimated  arrival  rates,  and  to  implement  the 
optimal  packet  transmission  probability  distribution  Pi  and  P2. 
Let  the  time  scale  of  network  dynamics  be  denoted  as  fast  time 
scale  (denoted  by  Time )  and  the  time  scale  of  arrival  statistics 
estimation  be  denoted  as  slow  time  scale.  In  this  experiment, 
a  slow  time-scale  epoch  contains  a  window  of  2000  fast  time- 
scale  epochs. 

Fig.  9  shows  the  steady-state  position  of  the  network  in  the 
equilibrium  phase  diagram  at  different  slow  time- scale  epochs. 
The  phase  diagram  is  the  same  as  the  one  shown  in  Section  III. 
Also,  the  time-series  response  of  D\  and  D2  is  shown  in  Fig.  10 
along  with  the  corresponding  control  set  point  Pi . 

For  model  identification  purpose,  parameters  are  identified 
to  be  a  =  -0.0771,  b  =  1.0404,  and  c  =  -b/a  =  13.4942. 
Initially  ( Time  =  1),  the  network  is  at  state  1  (as  shown  in 
Fig.  9)  with  Ai  =  0.05,  A2  =  0.01,  and  Pi  =  P2  =  0.5.  Over 
slow  time-scale  epoch  1,  Ai  and  A2  are  estimated,  and  at  the 
beginning  of  slow  time-scale  epoch  2  ( Time  =  2000),  the 
network  moves  to  the  optimal  location  with  Pi  «  0.8  (state  2  in 


Fig.  9).  As  a  consequence,  the  very  small  ripples  noticed  in  the 
Di  time  series  in  slow  time-scale  epoch  1  disappear.  (Note  that 
this  is  a  Phase  Type  II  control.)  Although  lower  Pi  would  have 
been  permissible,  it  reached  Pi  =  0.8  to  optimally  move  away 
from  the  phase  boundary.  At  the  beginning  of  slow  time- scale 
epoch  3  ( Time  =  4000),  the  packet  arrival  statistics  is  changed 
to  Ai  =  0.01  and  A2  =  0.05  (state  3  in  Fig.  9),  and  it  can  be 
observed  that,  after  some  initial  transience  response,  D2  settles 
to  a  nonzero  value  during  slow  time- scale  epoch  3  whereas  Di 
remains  at  zero.  Upon  estimation  of  Ai  and  A2  over  slow  time- 
scale  epoch  3,  the  controller  takes  the  network  to  the  optimal  lo¬ 
cation  with  Pi  ^0.26  (state  4  in  Fig.  9)  at  the  beginning  of  slow 
time-scale  epoch  4  ( Time  =  6000).  As  a  consequence,  D2 
response  dies  down,  keeping  the  Di  =  0  response.  The  network 
is  still  in  Phase  Type  II.  Finally,  at  the  beginning  of  slow  time- 
scale  epoch  5  ( Time  =  8000),  the  packet  arrival  statistics  is 
changed  to  Ai  =  0.15  and  A2  =  0.15  (state  5  in  Fig.  9);  both  Di 
and  D2  responses  settle  in  a  nonzero  value  over  the  slow  time- 
scale  epoch  5.  However,  the  steady-state  value  of  D 1  is  higher 
than  that  of  D2  as  the  packet  transmission  probability  distribu¬ 
tion  is  still  in  favor  of  packet  type  2  (with  Pi  =  0.26  and  P2  = 
0.74).  From  the  control  point  of  view,  the  network  is  in  Phase 
Type  I,  and  the  controller  drives  the  network  to  the  optimal  lo¬ 
cation  with  Pi  «  0.8  (state  6  in  Fig.  9)  at  the  beginning  of  slow 
time-scale  epoch  6  ( Time  =  10  000).  Higher  value  of  Pi  com¬ 
pared  to  that  of  P2  is  due  to  the  higher  priority  of  packet  type 
1  over  packet  type  2  (oq  >0^2).  Consequently,  D\  settles  to  a 
lower  value  compared  to  D2  over  the  slow  time- scale  epoch  6. 

V.  Summary,  Conclusion,  and  Future  Work 

This  paper  presents  statistical  mechanics-inspired  analysis 
of  critical  phenomena  in  square-grid  wired  communication 
networks  and  validates  the  pertinent  results  on  a  simulation  test 
bed.  A  comprehensive  finite-size  scaling  analysis  has  been  per¬ 
formed  for  a  specific  network  structure,  where  network  analogs 
of  intensive  thermodynamic  variables  (e.g.,  order  parameter, 
temperature,  and  pressure)  are  introduced  to  unambiguously  ex¬ 
plain  the  underlying  concepts.  Phase  diagrams  are  constructed 
for  networks  transmitting  heterogeneous  packet  types  for  the 
control  of  the  network  system.  Control  strategies  are  synthe¬ 
sized  to  achieve  robustness  in  the  fully  uncongested  phase 
and  to  mitigate  the  worst  case  packet  drop  rate  in  congested 
phases.  Closed-form  solutions  for  optimal  control  strategies 
have  been  formulated  for  the  square-grid  communication  net¬ 
work  structure.  The  analysis  framework,  reported  in  this  paper, 
can  be  potentially  extended  to  ad  hoc  wireless  sensor  networks 
for  boundary  surveillance  problems,  where  spatially  distributed 
autonomous  agents  with  (possibly  multimodal)  sensing  capa¬ 
bilities  collaboratively  monitor  a  given  environment.  In  general, 
different  sensors  sense  the  environment  in  a  distributed  manner 
and  need  to  communicate  with  other  types  of  sensors  for  in¬ 
formation  fusion  to  make  a  comprehensive  decision.  Therefore, 
complex  multihop  routing  of  (possibly  multipriority)  informa¬ 
tion  packets  in  a  sensor  network  bears  significant  relevance.  In 
addition,  such  a  simple  network  structure  may  serve  as  subsets 
of  large  complex  networks,  where  a  thorough  understanding  of 
small  subsets  is  necessary  for  the  control  of  the  entire  network. 
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Fig.  10.  Time  series  observation  of  packet  drop  rates  for  packet  type  1  (D i)  and  packet  type  2  (D2)  and  transmission  probability  for  packet  type  1  (control 
input)  Pi . 


In  this  paper,  the  packet  arrival  statistics  is  estimated  in  a  cen¬ 
tralized  fashion,  and  a  global  packet  transmission  probability 
distribution  is  broadcasted  for  optimal  control.  In  this  context, 
the  (statistical  mechanics-inspired)  ensemble  approach  may 
allow  the  control  algorithm  to  be  executed  over  a  network  in  a 
communication-constrained  environment  via  sparse  statistical 
sampling  for  the  estimation  of  network  performance  parame¬ 
ters.  It  will  be  more  useful  if  each  node  takes  its  own  decision 
based  on  the  local  (its  own  or  of  its  neighborhood)  estimation  of 
packet  arrival  statistics  and  yet  the  global  objective  is  satisfied. 
Due  to  the  ensemble  approach  of  the  current  study,  such  a 
distributed  control  problem  is  not  intractable  as  well  and  will  be 
an  important  topic  of  future  research.  This  approach  will  also 
enable  the  routing  strategy  to  handle  different  packet  arrival 
statistics  in  different  locations  of  the  network. 

This  initial  study  develops  a  priority-based  queuing  mecha¬ 
nism  and  does  not  include  the  aspect  of  scheduling  mechanism 
for  temporal  management  of  packet  transmission.  Therefore, 
from  an  application  perspective,  the  dynamic  aspects  of  the 
network  behavior  need  to  be  addressed  in  future,  particularly  to 
handle  distributed  implementation  of  the  algorithm  and  varying 
packet  arrival  statistics. 

Apart  from  the  issue  of  distributed  realization,  the  following 
topics  are  recommended  for  future  research  from  the  per¬ 
spectives  of  stability,  performance  analysis,  and  decision  and 
control: 

1)  validation  of  the  theoretical  results  in  more  complex  and 
realistic  network  scenarios  with  various  features  (e.g., 
layered  architectures  and  industry- standard  protocols); 

2)  investigation  of  the  effects  of  network  topology  (e.g. 
rectangular  instead  of  square  grid),  packet  arrival  statis¬ 
tics,  and  distributions  of  source  and  destination  nodes  on 
performance  analysis  and  decision  and  control; 

3)  dynamic  analysis  for  convergence  and  stability  of  the 
network  system  as  an  augmentation  of  the  equilibrium 
behavior  analysis,  reported  in  this  paper; 


4)  handling  of  distributed  implementation  of  the  algorithms 
and  varying  packet  arrival  statistics. 

Appendix  A 

Min-Max  Optimization 

Let  gi,i  =  1,  2, . . . ,  K,  be  continuous  and  strictly  monoton- 
ically  decreasing  functions  defined  as 

9i :  M  -  [grn,grx] 

where  0  <  x\  <  1,  <^(0)  =  g™ax ,  and  gi(xf)  =  g™171. 

The  optimization  task  is  to  minimize  the  cost  functional 

C(x  i,x2,  •  •  .,xK)  =  max  gi(xi) 

under  the  equality  constraint  Yhi  Xi  =  1  in  the  domain 
[x\X2  •  •  •  xK\  E  [0,  x\\  X  [0,  X$\  X  ...  X  [0,  xcK\.  It  is  noted 
that  the  feasible  region  is  nonempty  under  the  condition 
1.  In  this  context,  the  following  two  cases  are 
considered. 

Case  1)  Let  there  exist  a  solution  [x\  x\. . .  x*K\  (in  the  fea¬ 
sible  region  and  satisfying  the  equality  constraint) 
such  that  gi(x{)  =  52(^2)  =  •  •  •  =  9k(x*k)  =  9*> 
which  implies  that  the  cost  C{x\ ,  x2 , . . . ,  xk)  =  9 * • 
Now  let  a  nontrivial  perturbation  to  another 
point  in  the  feasible  region  be  [x\  x 2  ...  %]. 
Due  to  the  constraints  JA  x\  =  JA  Xi  =  1,  there 
must  exist  at  least  one  m  such  that  x *m  >  xm. 
Consequently,  <  9m{xm)  implying  that 

C(x\,x . . .  ,x*K)  <  C(x i,x2, . . .  ,xk)-  This  es¬ 
tablishes  that  [x\x*2  . . .  x*K]  is  the  globally  optimal 
solution. 

Case  2)  If  a  solution  of  the  form  [x\  x\  ...  x*K\  such 
that  gi(x\)  =  g2(x*2)  =  . . .  =  gK(x*K)  =  g *  does 
not  exist,  then  it  implies  that  the  solution  lies  on 
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the  edge  of  the  feasible  region.  A  constraint  is 
activated  by  choosing  an  index  m  with  the  smallest 
g™x  and  constraining  xm  =  0.  The  solution  of  the 
remaining  problem  is  similar  to  the  original  problem 
except  that  it  has  one  fewer  function.  This  solution 
is  locally  optimal,  and  the  process  is  iterated  recur¬ 
sively  to  arrive  at  the  final  solution.  In  other  words, 
the  scheme  is  to  sequentially  activate  constraints 
Xi  =  0  until  a  solution  is  found.  Consequently,  the 
maximum  number  of  iterations  that  may  be  required 
for  convergence  is  equal  to  K. 

Appendix  B 

Max-Min  Optimization 

Let  hi,  i  =  1,  2, . . . ,  K,  be  continuous  and  strictly  monoton- 
ically  increasing  homeomorphic  functions  defined  as 

hi  :[xli\^[0,h?ax] 

where  0  <  x\  <  1,  hi{x ?)  =  0,  and  hi(  1)  =  h™ax . 

The  optimization  task  is  to  maximize  the  cost  functional 

0[x i,  #2,  •  •  • ,  xk)  =  mmhi(xi) 

under  the  equality  constraint  JA  x^  =  1  in  the  domain 
[xi  X2  ...  xk]  £  [xi,  1]  x  [#2, 1]  x  . . .  x  [xcK,  1].  It  is  noted 
that  the  feasible  region  is  nonempty  under  the  condition 
J2ixi  <  !• 

Let  there  exist  a  solution  [x\  x\  ...  x*K]  (in  the  fea¬ 
sible  region  and  satisfying  equality  constraint)  such  that 
h\{x\)  =  ^2(^2)  =  •  •  •  =  =  h*,  which  implies  that 

0(X!,X2,  •  •  •  ,xk)  =  h*.  Let  a  nontrivial  perturbation  to  an¬ 
other  point  in  the  feasible  region  be  [x\  X2  ...  xk]-  Due  to 
the  constraints  JA  x\  =  JA  =  1,  there  must  exist  at  least 
one  m  such  that  x*m  >  xm.  Consequently,  hm(x:^n)  > 
implying  that  0(x{,  x\, . . . ,  x*K)  >  0(x\,  X2, . . . ,  xk)-  This 
establishes  that  [x\x*2  . . .  x*K\  is  the  globally  optimal  solution. 
The  subsequent  part  shows  the  existence  of  such  a  solution. 

Let  yer\f=i[0,hrx}  =  [y  mini  Vmax\,  where  ymin  =  0 
and  i/max  =  min^  h™ax .  Let  a  function  J  be  defined  on 
[0,  min*  h™0,00]  as 

K 

J(y)  =  ( y )  y  £  \y  mini  y  max]' 

In  the  above  equation,  the  inverse  exists  due  to  the  continuity 
and  strict  monotonicity  property  of  hi.  Under  the  assump¬ 
tion  of  homeomorphism,  h^1 ,  i  =  1,  2, . . . ,  K,  are  also  strictly 
monotonically  increasing  functions. 

The  evaluation  of  J  at  ymin  and  ymax  yields 

K 

JiVmin )  =  =  J2Xi  <  1 

i=  1  i 

K 

J(ymaX)  =  E  K1  (min  hr*)  >  ft-1  ( hZax )  =  1 


where  m  =  arg  min^  h™ax .  Therefore,  the  inequalities 
J  (ymin)  <  1  <  J  (ymax)  imply  that  there  exists  an 
h*  £  [ymin i  V max]  such  that  J(h*)  =  1.  The  optimal  solution 
is  obtained  as  x*  =  h^fo*),  which  satisfies  the  necessary 
condition  h\(x\)  =  ^2(^2)  =  •  •  •  =  hK{x*K)  =  h*  under  the 
constraint  JA  x\  =  \. 
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